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Nuclear dynamics in dense hydrogen, which is determined by the key physics of large-angle scattering or 
many-body collisions between particles, is crucial for the dynamics of planet's evolution and 
hydrodynamical processes in inertial confinement confusion. Here, using improved ab initio path-integral 
molecular dynamics simulations, we investigated the nuclear quantum dynamics regarding transport 
behaviors of dense hydrogen up to the temperatures of 1 eV. With the inclusion of nuclear quantum effects 
(NQEs), the ionic diffusions are largely higher than the classical treatment by the magnitude from 20% to 
146% as the temperature is decreased from 1 eV to 0.3 eV at 10 g/cm 3 , meanwhile, electrical and thermal 
conductivities are significantly lowered. In particular, the ionic diffusion is found much larger than that 
without NQEs even when both the ionic distributions are the same at 1 eV. The significant quantum 
derealization of ions introduces remarkably different scattering cross section between protons compared 
with classical particle treatments, which explains the large difference of transport properties induced by 
NQEs. The Stokes-Einstein relation, Wiedemann-Franz law, and isotope effects are re-examined, showing 
different behaviors in nuclear quantum dynamics. 

Nuclear quantum nature, especially for protons, manifests its significant effects on the ionic structures, 
phase diagram and ionic dynamics 15 . A lot of researches focus on the nuclear quantum effects (NQEs) 
such as quantum fluctuations and zero-point motion on the structures of dense hydrogen 610 , revealing 
many intriguing physics. However, the nuclear quantum dynamics such as ionic quantum transport is rarely 
studied because of the computational and the theoretical challenges. Only at low temperatures, the nuclear 
quantum dynamics was investigated. For example, large deviation for ionic diffusion and thermal conductivity 
of para-hydrogen is found below 32 K induced by NQEs compared with the classical-particle treatment 11 . At high 
temperatures, in particular in warm dense regime (i.e., above 1000 K), ions are usually considered as distinguish- 
able classical particles by assuming the quantum fluctuations or quantum collisions are negligible. However, even 
for the collision energies up to 10 eV, it has been shown that elastic collision cross sections for protons and 
hydrogen atoms can be very different in the quantum and classical theory 12 , indicating the possible failure of 
classical treatment for nuclear motion in warm dense regime. 

Dense hydrogen plays a crucial role in understanding the material behaviors under extreme conditions 13,14 and 
has important applications in energy sources 15 and astrophysics 1617 . The associated ionic and electronic transport 
properties are the key to tracing the dynamics of capsule implosion of inertial confinement fusion (ICF) 15 , 
modeling plasma processes 18 and the interior structure and evolution of giant planets and exoplanets 19-21 . 
More exact transport properties such as thermal conductivities and electron-ion relaxation rates are required 
because they can remarkably alter the dynamics of planet's evolution 22 and hydrodynamical processes in ICF 23 . 
The core pressures of some newly discovered exoplanets are estimated even up to 19 Gbar 21 , while the temper- 
ature may be several thousand Kelvins. Thanks to the recent progress in experimental techniques, one could get 
access to this previously inaccessible regime of phase diagram of hydrogen through a laser-induced shock wave 
loading of precompressed samples 24,25 . However, theoretically understanding the nature of hydrogen under the 
ultra-high pressures is still rare and a great challenge. Density functional theory based ab initio molecular 
dynamics (MD) has been widely used to investigate transport properties such as ionic diffusion, electrical 
conductivity of dense matter at high temperatures 26 28 . In viewpoint of collision physics, both the ionic and 
electronic transport properties depend strongly on the scattering cross sections between particles. In usual ab 
initio MD simulations, the corresponding scattering cross sections are statistically obtained from classical particle 
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dynamics simulations, in which the quantum dynamics such as tun- 
nelling is absent. The role of quantum dynamics in dense hydrogen is 
of great interest and expected to give a new perspective on NQEs. 

Path-integral molecular dynamics (PIMD) with imaginary- time 
path discretization has provided a proper theoretical description of 
the nuclear quantum fluctuation and derealization 29,30 . Further, the 
adiabatic centroid PIMD is proposed to obtain the real-time 
dynamics of nuclei in the quasiclassical approximation 31 . Its accuracy 
has been extensively studied for the bound systems 32-33 . However, for 
the continuum states, for example quantum tunnelling through a 
potential barrier, whether the centroid PIMD is accurate is already 
an open question. Since the ionic dynamics is strongly dependent on 
the processes of both bound and continuum states, it is urgent to 
understand this issue, looking forward to reach progress. 

In this work, we further developed the adiabatic centroid PIMD 
approach and the accuracy was demonstrated from quantum tun- 
nelling of a particle through the potential barrier. Within this theor- 
etical framework, the nuclear quantum dynamics of dense hydrogen 
was investigated, showing more interesting intrinsic physics. Both 
the ionic and electronic transport properties were greatly affected by 
NQEs, which enables us to understand the properties of dense matter 
more profoundly. 
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Figure 1 | Temperature dependence of the self-diffusion coefficient (a) 
and shear viscosity (b) at 10 g/cm 3 obtained from centroid PIMD and 
MD simulations. The viscosities obtained from the Stokes-Einstein (SE) 
relation are also presented, (c) Examination of SE relation with respect to 
temperature. The effective atomic diameter in SE relation is defined by the 
position of the first peak of radial distribution function. 



Results 

Structure and ionic transport properties. For dense hydrogen, five 
state points along the 10 g/cm 3 isochore from 0.1 eV to 1 eV, and 
three state points along the 1 eV isotherm from 10 g/cm 3 to 100 g/ 
cm 3 were studied. The main motivation for choice of these 
parameters is that in ICF a solid cryogenic shell of DT fuel in 
capsule needs to pass through a low temperature and high density 
stage before achieving the final state of the high density and 
temperature 34 . The transport properties at this stage are essential to 
understand the whole implosion process. Besides, the cores of some 
exoplanets are likely to be under the conditions studied here 13 . The 
structural properties were calculated using the primitive scheme 29 of 
PIMD, and the real-time quantum dynamics of nuclei was obtained 
through the improved ab initio centroid PIMD (see Methods). 
Exchange effects between protons are neglected in this study 35 . 
Firstly, we calculated the self-diffusion coefficient and the shear 
viscosity at 10 g/cm 3 with temperatures from 0.1 eV to 1 eV and 
the results are shown in Fig. la and lb. When taking into account 
the quantum collision properly, the ionic diffusion is much higher 
than that with classical treatment for ions. We can see that the self- 
diffusion coefficient at 0.1 eV from the centroid PIMD simulations is 
two orders of magnitude larger than the MD value. This large 
increment is likely to arise from the phase transition from solid to 
liquid with the inclusion of NQEs, which is confirmed directly by 
ionic trajectory analyses (see the Supplementary Information for 
simulated movies). The difference introduced by NQEs is also 
clearly demonstrated by radial distribution functions (RDFs) of 
hydrogen nuclei in Fig. 2a. Note that the RDF of hydrogen nuclei 
with MD method exhibits long-range ordered solid-like character at 
the low temperature of 0.1 eV, which is not visible in RDF obtained 
from the PIMD calculations. In the quantum-mechanical 
perspective, it can be understood from that the protons tunnel 
through the energy barriers and move away from their equilibrium 
positions in crystal structure at 0.1 eV, as appeared in hydrogen- 
bonded systems at low temperatures 5 . Our results accord with the 
conclusions in Ref. 36, i.e., the melting temperature of hydrogen is 
lowered by NQEs. In addition, we note that although hydrogen is in 
solid state at 0.1 eV in classical simulation, the self- diffusion 
coefficient of protons does not vanish, which arises from the 
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Figure 2 | Radial distribution functions of hydrogen nuclei from PIMD 
(solid lines) and MD (dashed lines) simulations with different 
temperatures at 10 g/cm 3 (a), and with different densities at 1 eV (b). The 

isotopes effects of deuterium (D) (long-dashed lines) and tritium (T) (dot- 
dashed lines) are also presented. The inset is the first peak of radial 
distribution functions at 100 g/cm 3 and 1 eV. 
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presence of atomic diffusion in solid hydrogen when the temperature 
is close to melting point 37 . 

The self-diffusion coefficients obtained via centroid PIMD are 
substantially larger than the classical particle value over the whole 
temperature range of 0.1-1 eV, as shown in Fig. la. The difference is 
from 146% at 0.3 eV, 61% at 0.5 eV to 40% at 0.7 eV. There is still a 
large difference of 20% even at the high temperature of 1 eV. It is not 
open-and-shut because the RDFs from the MD and PIMD calcula- 
tions are very close at 0.3 eV and almost identical when the temper- 
ature is up to 1 eV (see Fig. 2a). We note that the large-angle 
scattering between ions is dominant at such high densities consid- 
ered here. The pronounced quantum nuclear wave scattering cal- 
culation leads to a smaller large-angle scattering cross section 
between ions than the classical particle scattering calculation, and 
thus increases the mean free path of ions. Therefore, the ions diffuse 
more easily and the classical particle treatment of protons substan- 
tially underestimates the ionic diffusion. 

The shear viscosities v\ calculated from the Green- Kubo relation 27 
at the density of 10 g/cm 3 are shown in Fig. lb. Corresponding to the 
diffusions D, the viscosities from centroid PIMD are smaller than 
those from MD, but they become closer at high temperatures. In 
addition, the viscosities are also deduced from the Stokes-Einstein 
(SE) relation with slip boundary condition 38 . As discussed in Ref. 39, 
the SE relation, which is exact for the Brownian particles, is not valid 
in strongly coupled regime. It is shown in Fig. lc that the values of 
Dr\alk B T (a is the effective atomic diameter and defined by the posi- 
tion of the first peak of RDF) from both MD and PIMD calculations 
are much higher than the SE relation value 1/271. At low tempera- 
tures, pronounced NQEs induce more deviation from SE relation, 
indicating that the correlated time of particles in PIMD simulations 
(non-Markov dynamics) is longer than that of classical treatment. 

The largely reduced and broadened first peak of RDF from PIMD 
calculations at 0.1 eV and 10 g/cm 3 indicates the significant nuclear 
quantum derealization, which becomes weaker with increasing tem- 
perature and cannot be observed at 1 eV finally. Besides, RDFs show 
much pronounced nuclear quantum character with increasing den- 
sity (see Fig. 2b). When the density is increased up to 100 g/cm 3 , 
there is a distinct difference of RDFs between the PIMD and MD 
calculations due to NQEs even at the high temperature of 1 eV. 
Meanwhile, the isotope substitutions provide a clear demonstration 
of isotope effects on nuclear spacial distribution of condensed hydro- 
gen. Since the tritium nuclei is the heaviest one and has the shortest 
thermal de Broglie wavelength, RDFs of the tritons are more struc- 
tured than those of the other two isotopes. 

Electronic transport properties. Of our particular interest is 
whether the NQEs have evident effects on the electronic transport 
properties as profound as nuclei exhibits. To shed light on this issue, 
we calculated the electrical conductivity and optical absorption 
coefficient of dense hydrogen. In order to avoid the artificial drop 
at low frequency due to the finite number of atoms in the simulations, 
we employed the Drude formula to estimate the dc conductivity at a 
zero frequency limit 26 . It can be clearly seen from Fig. 3 that with the 
inclusion of NQEs, the low-frequency dependent electrical 
conductivity ff(co) exhibits a different trend compared with the 
results from the MD simulations, which directly leads to the 
different dc conductivity cr dc (w — » 0) except for the state point of 
10 g/cm 3 and 1 eV. The dc conductivities of dense liquid hydrogen 
from the PIMD calculations are largely suppressed compared with 
the MD value by 25% at 10 g/cm 3 , 0.3 eV and 18% at 100 g/cm 3 , 
1 eV. A visible larger optical absorption for hydrogen of 10 g/cm 3 at 
0.3 eV and 100 g/cm 3 at 1 eV is displayed in Fig. 3. It should be noted 
that the differences at 10 g/cm 3 and 0.1 eV originate from the 
different structures (solid in classical simulation and liquid in 
quantum simulation mentioned above). It is well known that 
NQEs will introduce more disorder for the systems, lowering the 
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Figure 3 | Electrical conductivity and optical absorption coefficient from 
the PIMD (solie lines) and MD (dashed lines) calculations for four state 
points, i.e., (10 g/cm 3 , 0.1 eV) (upper panel), (10 g/cm 3 , 0.3 eV) (second 
panel), (10 g/cm 3 , 1 eV) (third panel), (100 g/cm 3 , 1 eV) (lower panel). 

The error bars denote the standard deviation of averaging 10 atomic 
configurations along the molecular dynamics trajectory. Note that the 
structure from MD simulation at the state point of (10 g/cm 3 , 0.1 eV) is 
solid. 

melting point and critical point 9 . In fact, NQEs introduce ionic 
derealization and closer ionic distances, as shown in Fig. 2, which 
produce more localized electrons around the nuclei, evidenced by the 
distribution of charge density and electron localization function (see 
Fig. 4). The stronger localized electrons might enhance the optical 
oscillator strength and optical absorption cross sections of relative 
high frequencies. In addition, NQEs cause nonuniform broadening 
of energy bands because of the anharmonic effects, resulting 
significant difference in the electronic density of states (DOS). In 
particular, the DOS around Fermi level are decreased by NQEs 
(See the Supplementary Information). 

In Fig. 5, it is shown that the electronic thermal conductivities with 
the inclusion of NQEs are smaller than the value obtained from MD 
simulations as the electrical conductivity exhibits. We also computed 
the Lorenz number which is the ratio between electrical and thermal 
conductivity divided by the temperature 40 . According to the 
Wiedemann-Franz law, the Lorenz number is constant and reaches 
2.44 X 10~ 8 flWK -2 for the degenerate and coupled plasma 41 . 
Compared to the good agreement with Wiedemann-Franz law for 
the MD results, the Lorenz number with the inclusion of NQEs 
deviates slightly from the ideal value at low temperatures, where 
the ionic structures are significantly different. 

Discussion 

We have investigated the nuclear quantum dynamics in transport 
behaviors of dense hydrogen using the improved centroid PIMD 
simulations. The complex transport behavior shows that even when 
the NQEs affect little on the static structures (10 g/cm 3 and 1 eV),the 
ionic diffusions are still largely higher due to the lower collision cross 
sections with the inclusion of NQEs. In addition, electrical and ther- 
mal conductivities have a pronounced decrease in quantum simula- 
tions as well. The SE relation is not valid in strongly coupled regime, 
and the Lorenz number deviates from the degenerate limit at low 
temperatures. We have shown that the quantum nuclear character 
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Figure 4 | The cut-plane of 3D distribution of charge density (a, b) and electron localization function (c, d) of randomly selected configurations from 
MD (a, c) and PIMD (b, d) simulations at 10 g/cm 3 and 0.3 eV. In PIMD calculations, the configuration is randomly selected from imaginary-time slices 
in PIMD time steps. Note that with the inclusion of NQEs, the charge density and electron localization function have the larger maxima and smaller 
minima compared to the classical simulations. It indicates that electrons surrounding ions distribute more localized with quantum nuclei than classical 
particle treatment. All the other randomly selected configurations exhibit the same behavior as the example presented here. 

induces complex behaviors for both the ionic and electronic trans- 
port of dense hydrogen. 

Methods 

Centroid path-integral molecular dynamics for quantum tunneling. In the 

primitive PIMD formalism 29 , the N-particle quantum system is isomorphic to a 
classical system consisting of N ring polymers of length P with harmonic 
intrapolymeric forces. Thus the static properties of the quantum system can be 
obtained through MD simulations of its isomorphic classical system. In the adiabatic 
centroid PIMD scheme 31 , the normal mode transformation diagonalizes the 
harmonic bead coupling. The centroid and non-centroid modes are mass-scale 
separated. The nuclear quantum dynamics is obtained in the quasiclassical 
approximation through driving the centroid move in real-time in the centroid 
effective potential generated by all non-centroid modes. Taking for example the case 
of one kind of particle in the system, the coefficient of stiffness of harmonic 

interaction is cop — \/ r pksT— — VpEq, where P is the number of beads in the ring 
g 

polymer, k B is the Boltzmann constant, T is the target temperature of the simulated 
system, g is the total degrees of freedom of ions, £ 0 is the average kinetic energy of each 
ion determined by the target temperature T. In this scheme, cop keeps constant for a 
given temperature T regardless of the circumstance of the ions. 

Here we propose a new scheme in order to make the size of ring polymer is more 
close to ). as the instantaneous kinetic energy of ions is decreased. For each ion, cop 
depends linearly on its instantaneous kinetic energy E(t) during simulations, 
2N r 

Gjp(f) = — v P(aE(t) + bEo), where t is the simulation time step, two parameters a 




0.4 0.6 
T(eV) 

Figure 5 | Electrical conductivity, thermal conductivity and Lorenz 
number of dense hydrogen as a function of temperature at 10 g/cm 3 . For 

Lorenz number, the dot-dashed line denotes the value of the degenerate 
limit 2.44 X 10~ 8 QWKT 2 . 
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and b are determined by two boundary conditions. The first one is when E( t) = £ 0 , cop 

2N r- 

equals to the value obtained from the traditional formulation, i.e., — v PEq, resulting 

in a + b = I; the second one is co P reaches its minimum when E(t) — 0, and this 
minimum makes the radius of gyration of ring polymer up to half of the de Broglie 
wavelength of the corresponding ion. In PIMD simulations, the initial ring polymer is 
generated according to the Levy-Brown procedure 42 , and the corresponding radius of 
gyration of ring polymer R g oc 1 / vT. Since cop oc T, we have R* oc 1 /cop . Then a and b 
can be obtained according to the second boundary condition, that is, b = {IR^Il) 2 , a 
= 1 — {IRgolA) 2 , where Rgo is the radius of gyration of ring polymer when E(t) = E 0 , k 
is the de Broglie wavelength of the corresponding ion. In this way, the size of ring 
polymer is more close to a as the instantaneous kinetic energy of ions is decreased. 
The importance of this new scheme can be illustrated by considering the dynamics of 
a hydrogen atom tunnels through a Gaussian potential barrier using the improved 
centroid PIMD approach. The dependence of the transmission probability on the 
ratio of the incident energy E to the height of barrier U 0 is in good agreement with the 
WKB semi- classical approximation 43 and the exact results through solving the 
Schrodinger equation when a is less than the half-height width of the potential barrier 
(see Supplementary Information). Therefore, this new scheme is more suitable for 
simulating the nuclear quantum dynamics such as ionic diffusions, since they are 
dependent on the same physics of scattering or collision. 

Electronic structure calculations and MD simulations. All the PIMD and MD 
simulations were performed using modified Quantum-ESPRESSO package 44 . A 
periodic supercell including 250-432 atoms was employed according to different 
densities, which can ensure convergence of both ionic and electronic properties with 
good accuracy {see Supplementary Information for convergence tests). The 
Coulombic pseudopotential 45 and the generalized- gradient approximation 46 for 
exchange- correlation potential were adopted. The kinetic energy cutoff was set from 
200 to 400 Ry according to different densities. Electronic distribution was subject to 
Fermi-Dirac statistics 47 . The Brillouin zone was sampled using the F point in MD 
simulations. Denser k-point grids were used to calculate electronic properties. 

Within the framework of PIMD, the structural and thermodynamical properties 
were calculated using the primitive scheme 29 , while the real-time quantum dynamics 
of nuclei was obtained through the ab initio improved centroid PIMD. The self- 
diffusion coefficient was calculated from the slope of the centroid mean square 
displacement. The shear viscosity was obtained from the autocorrelation function of 
the off-diagonal components of the stress tensor 27 . Langevin thermostat 48 was 
employed to overcome the nonergodic problem, which not only produces a canonical 
ensemble and compensates the calculated errors 45 , but also gives us an efficient 
unified description from cold condensed matter to hot dense regime 45 ' 49 . Langevin 
thermostat was applied only to each noncentroid degree of freedom because ther- 
mostats would disturb the dynamical properties of the centroids. The Trotter number 
was set to 16 after a convergence test. Time step of 2-3 a.u. was used in our ab initio 
PIMD calculations with 40000 steps after thermalization, while the smaller time steps 
of 0.2-0.3 a.u. with 5 X 10 5 steps and a centroid adiabaticity parameter of 20 were 
adopted in the centroid PIMD simulations in order to decouple the centroid mode 
from other normal modes. 

The electrical and thermal conductivity and optical absorbtion coefficient were 
calculated via the Kubo-Greenwood formula 5051 with ABINIT package 52 . The elec- 
trical conductivity and optical absorption coefficient were calculated with a 0.03 a.u. 
width of Gaussian smearing and k-point sets of 6 X 6 X 6 and 8 X 8 X 8 for 10 g/cm 3 
and 100 g/cm 3 , respectively. Here, 10 atomic configurations along the trajectories 
were employed. In PIMD simulations, the electronic transport properties were 
obtained by averaging over both PIMD time steps and imaginary time slices. 
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